Screening in multilayer graphene 
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In this article we study the static polarization in ABC-stacked multilayer graphene. Since the 
density of states diverges for these systems if the number of layers exceeds three, screening effects are 
expected to be important. In the random phase approximation, screening can be included through 
the polarization. We derive an analytical integral expression for the polarization in both the full- 
band model and an effective two-band model. Numerical evaluation of these integrals are very time 
consuming in the full-band model. Hence, for ABC-stacked trilayer graphene, we use the two-band 
model to calculate the low momentum part of the polarization. The high momentum part is linear 
and is determined by calculating two points, such that we can determine the slope. The numerical 
results for the polarization of trilayer graphene are used to sketch the screened potential. 

PACS numbers: ... 



I. INTRODUCTION 

Stacking several layers of graphene on top of each other 
does not immediately lead to graphite. As long as the 
number of layers is small enough, the two-dimensional 
nature of the system is preserved, i.e. the (quasi) mo- 
mentum of the particles is oriented within the plane. The 
properties of these systems depend heavily on the way the 
layers are stacked and typically differ considerably from 
both monolayer graphene and graphite. 

There are two natural ways to stack graphene layers, 
namely AB or Bernal stacking and ABC (rhombohedral) 
stacking. In Bernal stacked multilayer graphene, the odd 
layers all have the same orientation, and so do the even 
layers. The orientation of the even layers is such that the 
B sublattice sites are opposite to the A sublattice sites of 
the layers directly beneath and above it. The A sublat- 
tice sites are located opposite to honeycomb centers. In 
rhombohedral stacked graphene, every layer is oriented 
such that the B sublattice is on top of the honeycomb 
centers of the layer beneath it and the A sublattice is on 
top of the B sublattice of the layer beneath it. This re- 
sults in a cyclic structure through different orientations. 
Hence, the layers i and i + 3 are exactly on top of each 
other. This lattice structure is shown in Fig. [I] 

Although a recent theoretical work investigates sys- 
tems in which the stacking of the layers is partly Bernal 
and partly rhombohedral, 1 so far most of the effort has 
been put into understanding either completely Bernal 
or completely rhombohedral stacked multilayer samples. 
These two systems behave very differently. In Bernal 
stacked multilayers, there are multiple low energy bands, 
i.e. quasi particles with different dispersions. When the 
number of layers is even (N = 2n), the n low energy con- 
duction bands are all parabolic (bilayer-like), but with 
different effective masses, while for an odd number of 
layers (N = 2n + 1) a linear band with the same slope 
as the energy band in monolayer graphene exists next to 
the n parabolic onesP The valence bands are related to 
the conduction bands by particle hole symmetry. On the 




Figure 1: (color online) Atomic structure of 
ABC-trilayer graphene. 



other hand, for ABC-stacked multilayers, the low energy 
physics takes place on the sublattice sites on the outer 
layers that do not have a direct neighbor in the next layer. 
As a result, it is possible to construct an effective 2x2 
Hamiltonian that is valid for energies E « t± ~ 0.3 
eVP From this effective model, it is easy to derive that 
the energy bands at small momenta and low energies dis- 
perse as E rsj k N , i.e. the bands become very flat when 
N increases. At the if -point, where the conduction and 
valence bands touch, the dispersion of the bands results 
into a diverging density of states when N > 3. This is 
in sharp contrast to Bernal stacked graphene, where the 
density of states never diverges at the Dirac point. 

The integer quantum Hall effect could be a way to 
identify the different stacking orders. This is because 
the Landau level spectrum is very different in the two 
systems W For example, in the trilayer case the Landau 
levels disperse with the magnetic field B as E ~ B 3 ^ 2 for 
rhombohedral stacking, while the linear and parabolic 
bands in Bernal stacked trilayers give rise to two sets of 
Landau levels. One set disperses as a graphene mono- 



layer, E m \ ~ \[B, while the other behaves as E\>\ ~ B, 
just as a graphene bilayer does. Hence, the Landau lev- 
els cross as a function of the magnetic field. 5 This fun- 
damental difference is true for any N > 3 multilayer 
system: When the stacking is Bernal the Landau lev- 
els cross, while for ABC-stacked systems they do not (as 
long as the high-energy bands are neglected) . Due to the 
low mobility of most multilayer graphene samples, much 
higher magnetic fields are required to observe the quan- 
tization of the Hall conductance. Nevertheless, for both 
Bernal and rhombohedral stacked trilayer graphene the 
integer quantum Hall effect is observed in experiment P^ 
In the Bernal stacked case, hexagonal boron nitride was 
used as a substrate,^ increasing the mobility by a fac- 
tor of 100. This technique may be used in the future to 
observe the quantum Hall effect in graphene multilayers 
with an even higher number of layers. 

Although the integer quantum Hall effect was observed 
promptly, it took much longer to confirm that the frac- 
tional quantum Hall effect exists in graphene. As a result, 
the importance of the Coulomb interaction in graphene 
has long been debated. Theoretical predictions of inter- 
action effects had been made, such as a ferromagnetic 
phase transition in monolayer, bilayer, and later also in 
trilayer grapheneP^ The observation of the fractional 
quantum Hall effect in 2009 confirmed that interactions 
do play a role in graphene physics. The groups of Eva An- 
drei and Philip Kim reported the quantum Hall plateau 
in suspended graphene with high mobility at filling fac- 
tor v = 1/3, using a two-terminal deviceP^^ However, 
a two-terminal setup cannot provide an unambiguous 
proof of the existence of the phenomenon. The issue has 
only been definitively settled after a four-terminal device 
was used to observe the v — 1/3 plateau in suspended 
graphene^ and several other plateaus at fractional fill- 
ing factors in graphene on hexagonal boron nitride. 16 The 
importance of interactions in graphene was later reiter- 
ated by other experiments, for example the renormaliza- 
tion of the Fermi velocity due to the Coulomb interaction 
in monolayers. 17 In addition, there is evidence that the 
fractional quantum Hall effect also occurs in suspended 
bilayer and trilayer samples P^ 

The Coulomb interaction is always present in systems 
that consist of many charged particles, like electrons. 
The importance of electromagnetic interactions depends 
on the properties of the system. The ratio of the Coulomb 
to kinetic energy r s is a measure of the influence that the 
Coulomb interaction has on the system. When the ki- 
netic energy dominates and r s is small, the system can 
be described as a Fermi liquid. When r s is large new 
phases can occur. For a two dimensional electron gas 
(2DEG), r s = m*e 2 /(e^ 2 v /7rn e i), where m* is the ef- 
fective mass of the electrons, e the electron charge, e 
the dielectric constant, and n e \ the density of electrons. 
Hence, for low electron densities the Coulomb interaction 
dominates and other phases, for example a Wigner crys- 

— 1/2 

tal, can form. The n el dependence is the result of the 
Coulomb interaction (V) ~ l/(r) ~ v^I and a quadratic 



kinetic energy (K) ~ k\ ~ n e \. When the dispersion is 
not quadratic, r s will become a different function of the 
electron density. 

In monolayer graphene, the charge carriers behave as 
massless relativistic particles. Therefore, the kinetic en- 
ergy scales with momentum or y / n^, instead of momen- 
tum squared or n e \ as it was the case in the 2DEG. Hence, 
the parameter r s depends only on material parameters 
and is independent of electron doping r s = e 2 /(efti^)- 
For graphene, e is the average dielectric constant of 
the material below and above the graphene layer, i.e. 
e = 1 for suspended graphene in vacuum and e = 2.5 
for graphene on a Si02 substrate. Thus, r s = 2.2 and 
r s = 0.8 for these two cases, respectively. Compared to 
a typical 2DEG, graphene is weakly interacting. How- 
ever, close to the charge neutrality point the density of 
states vanishes. As a consequence, there are not many 
electrons available for screening and the Coulomb inter- 
action is almost unscreened, thus remaining long ranged. 
Indeed, the Thomas-Fermi screening vector, which is the 
k — >• limit of the polarization, scales with the Fermi en- 
ergy and therefore vanishes if the system is close to half 
filling. 19 In the short- wavelength limit, where k is large, 
the polarization is linear and, in this regime, the effect of 
screening is a renormalization of the interaction strength. 

For bilayer graphene, the parameter r s scales as r s ~ 
1/v^ei- 19 Hence, close to half filling, where n e \ = 0, the 
interaction term should dominate the kinetic term. How- 
ever, not only is it very difficult to produce a charge 
neutral system, due to the formation of electron hole 
puddles,^! it is also no longer true that the Thomas- 
Fermi vector vanishes for n e \ = 0. The Thomas- Fermi 
vector is independent of the density of electrons in bilayer 
graphene. Therefore, screening is more profound in bi- 
layers than in monolayers of graphene. The polarization 
can be calculated analytically and from the polarization 
it is possible to construct the screened potential.^ Due to 
the two-dimensional nature of the system, the potential 
is not exponentially screened, but remains polynomial. 

For ABC-stacked multilayers with three or more lay- 
ers, the density of states diverges at the charge neutral- 
ity point. Since the Thomas- Fermi vector scales with the 
density of states, it is expected that screening is impor- 
tant for such systems. Nevertheless, not much is known 
about Coulomb interations in multilayer graphene sys- 
tems. 

The aim of this paper is to determine the polarization 
and the screened potential in rhombohedral stacked mul- 
tilayers. Firstly, two models are introduced in section 
HH In the full-band model the full 2N x 2N Hamiltonian 
is used, while in the two-band model an effective 2x2 
matrix is introduced. In section |III[ the polarization is 
calculated in the two-band model and it is shown that 
this approximation breaks down for large momenta. Al- 
though we analytically derive the formal integrals which 
have to be solved to calculate the polarization in N lay- 
ers of graphene, we solve the problem numerically only 
for ABC-stacked trilayer graphene, as an example. We 



also show results for the full-band model. The screened 
Coulomb potentials are derived for the two-band model 
in section IV and a realistic sketch of the screened po- 
tential is drawn in the full-band model. We discuss our 
results in section IVl 



II. THE MODEL 

A. Full-band Hamiltonian 

To describe an ABC-stacked multilayer of graphene 
with N layers, we use a nearest-neighbor tight-binding 
model. Hence, the electrons can tunnel to adjacent lat- 
tice sites within the same layer (with energy t = 3 eV) 
or to direct neighbors at a distance d = 3.4 A in other 
layers (with energy t± = 0.35 eV). The noninteracting 
tight-binding Hamiltonian in real space is given by 
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If c € {a, 6}, then c\, a (cj. j<7 ) creates (annihilates) an 
electron on lattice site i in layer I E {1,2, ...,iV} with 
spin ae {t,|}- 

Since the unit cell of this system contains 2N lattice sites, 
the reciprocal space representation of this Hamiltonian is 
a 2 N x 27V matrix. After expanding around the if-point, 
the low-energy Hamiltonian is cast into the form, 



H = fd 2 k^(k)n ^(k), 

n = 



{H ml B 
B T H m] B 



(1) 



B 






ke i(f) M 





t± 



^ t (k) = (a t 1 (k),6 t 1 (k),4(k),...,6j v (k)), 

where hvp = (3/2) at defines the Fermi velocity in mono- 
layer graphene. Furthermore, k is the norm of the two- 
dimensional momentum vector, k = |k|, and </>(k) = 
arctan (k y /k x ) is the angle of the momentum vector. In 
the following, we refer to Eq. [I] as the Noninteracting 
full-band Hamiltonian. 



B. Two-band Hamiltonian 

In an ABC-stacked multilayer, only the A sublattice in 
the bottom layer (layer 1) and the B sublattice in the top 
layer (layer N) do not have direct neighbors in an adja- 
cent layer. The electrons on sites with a neighbor in an 



opposite layer will dimerise and the energy bands associ- 
ated with these electrons will move away from the charge 
neutrality point. This results in two energy bands close 
to the charge neutrality point, while the other energy 
bands split away from zero by an energy ~ t±. Hence, for 
an ABC-stacked multilayer of graphene, the low-energy 
physics takes place on the A\ and the Bn sites. There- 
fore, it is possible to construct an effective low-energy 
model that takes only the two energy bands into account 
that are closest to the Dirac point P^ This low energy 
Hamiltonian is a 2 x 2 matrix and since it takes N intra- 
plane and TV — 1 inter-plane hoppings to go from the A\ 
to the Bn site, it has the form, 
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We will refer to the Hamiltonian in Eq. [2] as the Non- 
interacting two-band Hamiltonian. This Hamiltonian is 
valid for small momenta at which the energies are much 
smaller than t±. 



III. THE POLARIZATION BUBBLE 

The Feynman diagram of the polarization is shown in 
Fig. |2| In the random phase approximation, this bubble 
diagram can be used to compute the screened potential 
or the free energy in an infinite order expansion. The 
screened part of the potential can also be absorbed into 
the dielectric constant. By doing so, one can relate the 
polarization to the electromagnetic susceptibility x(cj, k), 
which is defined by e(cj, k) = l+47rx(^, k), and measures 
the tendency of the medium to adjust to an external elec- 
tromagnetic perturbation. Furthermore, the k — »■ limit 
of the polarization gives the Thomas- Fermi screening vec- 
tor. The dynamical part of the polarization is needed 
to describe plasmons, but those will not be treated in 
this paper. Here, we focus on the static polarization, for 
which lj = 0. 




Figure 2: The bubble diagram n^#- 



A. Two-band model 

In this section, we will calculate the polarization and 
subsequently the screened Coulomb interactions in the 
two-band model. We can neglect spin in this problem. 
The Hamiltonian is a matrix and hence, the bubble di- 
agram will have indices labeling lattice site. For conve- 
nience, we indicate the A\ sites by A and the Bn sites 
by B. The polarization is given by 



n 



lUa Rab 

12, 



/i'2 

x Gji(iQ n ,q), 



-icj m ,q + k) 



(3) 



where GV/(zco>,k) is the electron propagator between the 
lattice sites i and j, uo m and Q n are Matsubara frequen- 
cies, and T is temperature. In Fig. [2] the Feynman dia- 
gram of the Hab bubble is shown. 

The propagator, which is a 2 x 2 matrix in the two- 
band model, is given by G(zco>, k) = (iujt — %o B ) _1 , where 
1 is the identity matrix. A derivation of the propagator 
in the full-band model is given in the appendix, but can 
be applied in the two-band model as well. The diagonal- 
ization matrices of Hq B are defined as 
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Writing G(zcj m ,k) in this form allows us to perform the 
Matsubara sum in expression ([3]) for the polarization: 
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where is the angle between k and q and n(£ s ) is the 
occupation function of the energy band s. Note that 
cos (|(iV — ra)7r) G { — 1, 0, 1}. For fixed TV, we have de- 
termined the structure factor Ff s (k,q,0) analytically. In 
the two-band model, this factor depends on and on the 
ratio q/k. In the full band model we can only determine 
this factor numerically. Since it is (in that case) a func- 
tion of /c, <?, and #, it will slow down the calculations. 



where we have defined a = t±(hvp/t±) N . Let us denote 
£| = sak N , where s = ±. Then, 

A+ = diag(l,0), 
A~ =diag(0,l). 



Zero Fermi energy 

Let us first consider the half-filled system for which 
Ep = 0. In this case, the valence band is completely 
filled, while the conduction band is empty, 

»(#) = 0, n(4") = 1. 

Let us define x = q/k; then the expression Q for the 
polarization can be written as 
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To arrive at this expression, we made the denominator will yield a constant that can be determined numerically, 
real and filled in the expression for £|. It is now obvious 
that the polarization is real, as long as uj m = (static 
screening). For static screening, one may extract the k 
dependence of the polarization out of the integral, which 
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Figure 3: (color online) Plot of aIl AA (uj — 0, k) (blue dashed line), aU AB (oj = 0, k) (red dotted line), and 
aU® ot (uj = 0, k) (black solid line) in the two-band model for a different number of layers and kp = 0. 



For the ABC-trilayer, one finds that 
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We conclude that for an undoped trilayer, the polariza- 
tion goes as ~ 1/fc and therefore diverges as k — >• 0. 
This is expected, because the long wavelength limit of 
the bubble is proportional to the density of states at the 
Fermi energy and for an ABC- stacked trilayer the den- 
sity of states diverges at the Dirac point. This singular 



behavior is present for all ABC-multilayers with N > 3, 
as can be seen in Fig. [3] 



Nonzero Fermi energy 

Now, let us assume a nonzero Fermi energy Ep > 0. In 
this case, the conduction band is partially filled. Hence, 
the occupation functions become 

n(Z+) = Q(k F -k), n(4") = l, 

where Q(x) is the Heaviside theta function and kp is the 
Fermi momentum. If we write down the expression for 
the polarization, we recognize the expression for the half 
filled case, li^Ak) plus a correction 
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The term with F~-~ (1, x, 9) vanishes, because both oc- 
cupation functions are unity and therefore cancel each 
other out. From the expression above, we learn that 
only the extra terms with respect to the undoped case 
depend on the Fermi momentum. These correction terms 
contain two Heaviside theta functions. One of them, 
Q(kp/k — x), is nonzero only within a circle of ra- 
dius S = kp/k around the origin, while the other one, 
Q(kp/k — y/l + x 2 + 2xcos#), is nonzero within a dis- 
tance S from the point (1,0). Since 5^0 when k — » oo, 
these correction terms will go to zero for large momenta. 
We conclude that, within the two band model, the large 
momentum dependence of the polarization always equals 
that of the half- filled system. Note that for kp/k < 1/2, 
i.e. k > 2kp-, the two regions where the Heaviside func- 
tions are nonzero do not overlap. For k < 2k p these 
two regions do overlap and this explains the cusp in the 
polarization at k = 2kp. 

The long wavelength limit of the polarizations are 
given by 
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where gp> is the degeneracy factor due to spin and valley. 
This condition is indeed satisfied by the numerical results 
we have found for the trilayer case. In Fig.|4a|the compo- 
nents Hij(k) and H tot (k) are plotted for a nonzero Fermi 
momentum in the ABC-stacked trilayer. After scaling 
the axes as is done in the figure, the plot is universal, i.e. 
independent of Fermi momentum. 
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Figure 4: (a) The normalized polarization functions for 
ABC-stacked trilayer graphene in the two-band model: 
H-AA(k) (dashed line), U AB (k) (dotted line), and U tot (k) 
(solid line). This plot is universal as long as kp ^ 0. (b) 
Plot of aU tot (uj = 0, k = 0) as a function of k F . 



B. The full-band model 

For the full-band model, the expression for the polar- 
ization can be derived in a similar way as we did to obtain 
Eq. Q. The main difference is that all the matrices are 
now 27V x 27V and in general, cannot be calculated an- 
alytically anymore. A derivation for the propagator is 
shown in the appendix. We have that 
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where z, j now label the 27V different sublattice sites, s 
and s' label the 27V energy bands, and therefore A s are 
matrices with zero's everywhere, except for a 1 in the 



entry along the diagonal. U^ is the diagonalization 



matrix of Hamiltonian (HI). For the full-band model we 
define n tot (fe) = E??=i ^(fe). 

The computation of the static polarization (oo = 0) 
is very tedious if the number of layers is greater than 
three. We can argue how H tot (k) behaves in the short- 
and long- wavelength limit. The relation (pi still holds 
and therefore we can plot II (a; = 0, k = 0) as a function 
of Fermi momentum (Fig. 4b). As expected, in the small- 
k limit the total polarization of the full model agrees with 
the results we found in the two-band model. 
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Table I: Numerical values of U tot (k) in ABC-stacked 
trilayer graphene for kp = 0.017. 

Although the two band model describes the polariza- 
tion well for small momenta, the short-wavelength behav- 
ior differs dramatically. We found that in the two-band 
model II tot (/c — >• oo) ~ l/k N ~ 2 . This relation was in- 



duced by the k N dispersion. In the full band model, this 
dispersion relation is only valid for momenta close to the 
Dirac point. For larger momenta, the dispersion of the 
bands eventually becomes linear. If k is large, this lin- 
ear regime of the dispersion dominates the polarization 
integral and therefore II(fc — >• oo) ~ fc, as is the behavior 
for monolayer graphene. 23 Hence, the short- wavelength 
limit of the static bubble is linear. This linearity is inde- 
pendent of TV. The slope does depend on the number of 
layers, but not on the Fermi energy. For trilayer graphene 
we have confirmed the linear behavior for a fixed Fermi 
energy {kp = 0.017) through a straight forward numer- 
ical calculation. This very time consuming process re- 
sulted in two points that align perfectly with the origin 
(see Table pi), confirming that Ii tot {k — >• oo) — — 7/c, with 
7 w 0.18. In Fig. 1 " 



5a 



the two limiting regions of H tot (k) 
are shown. The low momentum behavior, computed in 
the two-band model, is universal after scaling. Since the 
slope of the linear part is fixed, it is not invariant after 
scaling and therefore depends on kp. Hence, the total 
polarization in the full-band model is no longer univer- 
sal. 

In this paper, we concentrated on the behavior of the 
polarization in the static limit (cj m — » 0). A generaliza- 
tion to include the dynamical part would be an interest- 
ing topic for further studies. 
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IV. THE SCREENED POTENTIALS 
A. Two-band model 

Now that we have obtained the polarization functions, 
it is possible to determine the screened Coulomb poten- 
tials. In the two band model we only have the potentials 
within the layers and the potential between the first and 
AT-th layer, since the low energy physics takes place on 
the A\ and Bn sublattice sites. It is difficult to estimate 
the results for the nearest-layer potentials, Vj il B 2 {k) for 
example, since we cannot say anything about the nearest- 
neighbor polarizations, e.g. TL Al & 2 . In fact, in the full- 
band model we need all the components of the polariza- 
tion matrix, if we want to calculate any screening poten- 
tial. Nevertheless, let us use the two band model for now 
and discuss its limitations later. 

The screened potentials are solutions of a Dyson equa- 
tion (see Fig. |6l). However, since the potentials carry 
two layer indices, the Dyson equation is a matrix equa- 
tion in this case. We focus now on the trilayer case. If 
one writes down the equations for the intralayer poten- 
tial VaxAiQz) — VB 3 B 3 (k) = VAA(k) and the interlayer 
potential V t/ \ 1 s 3 {k) = V AB {k) : one notices that they are 
coupled. By defining T as 

T(k) = {[V AA {k) - V AB {k)][TL AA {k) - U AB {k)} - 1} 

x {[V AA {k) + V AB (k)][U AA (k) + U AB (k)} - 1} , 



Figure 5: The solid lines sketch (a) the full-band 

polarization and (b) the full-band screened potential 

y tot (A:) for ABC-stacked trilayer graphene. The \ow-k 

regime is obtained in the two-band model and the 
high-A: regime is obtained through a direct numerical 

computation, (a) The full-band polarization for 

E F /a = 0.027 (k F = 0.3). (b) The full-band screened 

potential for E F /a = 0.001 (k F = 0.1). In the high-& 

regime the potential becomes a rescaled version of the 

unscreened potential V(k) = 1/fc. 



the solutions are given by 



V n (k) 

V 13 (k) 



V AA {k) - n AA (k)[V AA (k) 2 - V AB (k) 2 } 

T(k) 
V AB (k) - U AB (k)[V AB (k) 2 - V AA (k) 2 } 

T(k) 



WWW 
' j 



Figure 6: The Dyson equation, which is a matrix 
equation. 
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Figure 7: (color online) The screened potentials in the 
two-band model: Vaa (red dotted line), Vab (blue 

dashed line), and V tot (black solid line). The 

unscreened interaction V(k) = l/(fc) is given by the 

purple dotdashed line. The interlayer distance d is set 

equal to 1 and kp =0.1. 



where the bare interactions read 



V AA (k) 
V AB (k) = 



fc' 



-2kd 



k 



and d is the interlayer distance.' The screened poten- 
tials are shown in Fig. [7| For convenience we have also 
defined V tot (k) = l/[fc + n tot (fc)], which is the solution to 
the Dyson equations if we set all exponentials exp(— 2kd) 
equal to unity, which is a valid approximation when k 
is small. V tot (k) has the same features as the real so- 
lutions, but because of its much simpler form, it will be 
convenient to use this potential during more in depth 
discussions. 

Although this approximation describes well the screen- 
ing for small momenta, the large k behavior of the 
polarization in the full band model is different than 
the l/k N ~ 2 decay in the two-band model. Hence, the 
screened potentials in the two-band model will also be 
incorrect in the large- k limit. 



B. Full-band model 

The screened potentials in the full-band model are also 
solutions of the Dyson equation sketched in Fig. [6) which, 
for trilayer graphene, is now a 6 x 6 matrix equation. The 
solution is given by 

yscr = ( t _ VI iyi v 

It is easy to obtain the screened potentials numerically, 
once all components of the polarization are known. How- 
ever, as we have seen previously, the calculation of the 
components of n is numerically very time consuming in 



the full-band model. Nevertheless, we obtained a re- 
alistic sketch of n tot (/c) = T, 2 ij=i u ij( k )' As before, 
when we put the exponentials in the potentials to unity 
(equivalently put d = 0), the screened potential reduces 
to V tot (k) = l/[k + U tot (k)]. Therefore we can sketch 
V tot (k) as well, which is done in Fig. 5b 



V. DISCUSSIONS 

In this paper we have derived an expression for the 
polarization for ABC-stacked multilayer graphene. The 
calculations were performed within both the full-band 
model and the two-band model, an effective low-energy 
model in which the 27V x 27V matrices reduce in size to 
2x2. The advantage of the effective model is that it 
becomes easier to calculate the polarization numerically. 
The drawback is that the large-fc behavior of the polar- 
ization becomes flawed. Instead of the linear behavior 
n ~ —<yk at large fc, which is imposed by the linearity of 
the energy bands farther away from the Dirac points, the 
polarization drops off as l/k N ~ 2 in the two-band model. 
It is very time consuming to calculate the full-band po- 
larization numerically, such that the results of the two- 
band model are of great help to understand the behavior 
of n tot (fc) = Y^ N j=i u ij( k ) in the full-band model. One 
has to be aware that in the effective model only the A\ 
and Bn lattice sites are involved. Therefore, we can say 
nothing about the polarization functions between other 
lattice sites in this approximation. 

We have confirmed the linear behavior of n tot (fc) for 
the ABC trilayer by calculating two points in the full- 
band model at high momenta. These two points align 
with the origin, confirming the linear asymptote. Al- 
though we do not know the exact cross over between the 
low momentum and the linear regime, we can sketch the 
total polarization (see Fig. 5a). The curves obtained in 
the two-band model for the polarization are universal af- 
ter scaling, the cross over to the linear regime is not. This 
is because the slope of the linear regime is constant. Af- 
ter scaling however, the slope becomes — 7/n tot (0). Since 
|II tot (0) | is decreasing as a function of kp (see Fig. 4b), 
the linear asymptote will become steeper in the scaled 
plots when kp is increased. For example, the cross over 
to the linear regime is around 2fc^ for Ep/a = 0.03 
(kp = 0.3), while it is around Skp for Ep/a = 0.001 
(kp = 0.1). Hence, the polarization is no longer univer- 
sal after scaling in the full-band model. 

Comparing our results for the ABC trilayer with the 
full-band 21 and two-band 24 results for the bilayer, which 
are shown in Fig. [8j we see some similarities. As long 
as the Fermi energy is nonzero the total polarization is 
finite everywhere. The fc — >• limit is proportional to 
the density of states, just as in the bilayer case. This is 
in fact true for any TV-layer ABC-stacked sample. It is 
important to realize that, although the plots look simi- 
lar, the absolute values of n tot (fc — »> 0) are much higher 
for the ABC-stacked trilayer than for the bilayer. This 
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Figure 8: Rescaled polarization (II tot (fc)) for the 

graphene bilayer. Picture extracted from Ref. [2TJ /i is 

the Fermi energy and t equals t± in the notation used 

throughout this paper. 



is because the density of states is much higher. The nor- 
malized polarization has a peak as k increases and then 
a cross over to a linear regime takes place if k is further 
increased. This cross over is not captured by the two- 
band model. There are also differences in the behavior 
of bilayer and trilayer samples. The peak in the bilayer 
is located exactly at 2kp and has a discontinuity in the 
first derivative. For the trilayer the peak is smoother and 
the maximum is reached before 2k p. The discontinuity 
in the first derivative at 2k p remains, although it is less 
pronounced. Also for the bilayer the crossing to the linear 
regime is shifting to smaller momenta (measured in units 
of kp) when the Fermi momentum increases. The slope of 
the linear part is different for the trilayer compared with 
the bilayer (and the monolayer). This can be explained 
by the existence of more valence bands which are filled 
in the trilayer. In general, for an ABC-stacked iV-layer 
system, there are N — 1 filled valence bands (and also 
N — 1 empty conduction bands) further away from zero 
energy. These bands are expected to give a contribution 
when computing the full-band polarization. Hence, there 
is no reason that the slope of the linear asymptote should 
be independent of layer number N. We expect that the 
slope of the linear part of the total polarization increases 
further when N grows larger. The static polarization for 
monolayer graphene is constant up to k = 2k f and then 
becomes linear with a slope 7 m .i. = 1/16P3 Comparing 
this with the slope of the bilayer 7b. i. = 1/8 = 2/16, and 
the trilayer 7 = 0.18 ~ 3/16, a trend is seen. Although 
our results are not accurate enough, this is an indication 
that the slope of the linear part of the polarization scales 
linearly with the number of layers. 

The unscreened Coulomb potential V(k) = g/k with 
interaction strength g diverges as k — >• 0. If kp = 0, the 
polarization diverges as well for k —> 0. The screened 
potential will converge to zero in this limit, V tot (k — >• 
0) = 0. When kp =^ 0, the potential will be finite and 
nonzero everywhere. The screened potential has a lo- 



cal minimum and will converge to a renormalized version 
of the unscreened potential in the large- k limit. Due to 
the linear behavior of the polarization II tot (/c) = — jk 
at large momenta, the screened potential is just the un- 
screened one with a renormalized interaction strength in 
this regime. The new interaction strength has the form 
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1 + 7 



Since 7 is positive, the interaction strength will be re- 
duced. A sketch of the screened potential is shown in 
Fig." 



5b 



With the insights gained here, we conclude that the 
simplified model for the polarization proposed in Ref. [12] 
to include the effect of screening in the ferromagnetic 
phase transition in ABC-stacked trilayer can be further 
refined. Indeed, instead of using the slope for the bilayer, 
7 = 0.125, we have found that the precise value for trilay- 
ers should be 7 = 0.18. As a result, the critical couplings, 
and therefore also the critical doping levels, are reduced. 
The reduction of the interaction parameter g becomes 
larger when g increases. For g = 6 the extra reduction is 
around 28%. Since the phase boundary n CT [t(g) is quite 
fiat for g > 1, the critical curve of the screened case in 
the phase diagram will not change drastically. 

In sum, the results presented here allow one to deter- 
mine the effect of (static) screening in N layers of ABC- 
stacked graphene and should thus contribute to more 
accurate calculations of interaction phenomena in mul- 
tilayer graphene. The number of trilayer experiments 
is increasing over the last few years, hence theoretical 
work is needed as well. By connecting the universal low- 
momentum part and the known high-momentum linear 
regime of the polarization, our numerical data could be 
used to include screening in a more realistic manner in 
the description of these fascinating materials. 
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A. GREEN'S FUNCTION 

In order to derive the non-interacting Green's function 
for ABC- multilayer graphene, we can use the Feynman 
path integral formalism. Using the definitions in Eq. [T] 
the partition function is given by 



Z = 



y"d^t]d[# 



-s[ijj^ijj]/h 



where the action is given by 



S= f dr j d 2 k^(k) (h^+H (k)\ ^(k). 



In this expression, r denotes imaginary time. Next, de- 
fine C/k> such that D(k) = U^Ho(k)U\z is diagonal and 
unitary. /7k is called the diagonalizer of the Hamilto- 
nian. For the full-band model it can only be determined 
numerically. The diagonalizers induce a change of basis 
that is used later y?(k) = Uip(k). 



^k 



k,cr — ^ C k,cr,l> C k,cr,2> ' C k,cr,2iV-l ' C k,cr,2iV J J V>) 



where c^ a a (ck,o-,a) creates (annihilates) an electron 
with momentum k, spin <j in energy band a. The Hamil- 
tonian in Eq. (fil) can be rewritten as 



^o = 5>l,^(kW. 



(7) 



The Green's function is defined as (^''"(k)'^(k)), which 
is in fact equivalent to the inverse of the quadratic part 
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of the action. After performing a Fourier transforma- 
tion from imaginary time to Matsubara frequencies and 
neglecting spin, this results into 

G(iuj m , k) = [iu m l - H (k)] _1 , 

= [iLj m l-Uy,D(k)Ui\~\ 

= {U±[iu m ±- D(k)]uiy\ 

= ul[iuj m t-D(k)]- 1 U^ 

2N 

= E- J re^ f / k , 

s=1 IUJ ™ Sk 

where ££ are the diagonal entries of D(k), which are the 
eigenenergies of the Hamiltonian, and A s is a 2N x 27V 
matrix with zero's everywhere, except for a 1 in the 5 th 
entry along the diagonal. In other words, the sum over 
s is over the different energy bands. This Green's func- 
tion is a 27V x 27V matrix. The components correspond 
with sublattice index. Hence, the propagator between 
two different sites i and j is defined as 



Gij(iu) m ,k) 



2N 

E 



lUr, 



{ul^) i 
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